# Figure_2.R 

# Part of the replication archive for 
#
#   Bullock, John G. 2020. "Education and Attitudes toward Redistribution in
#   the United States." British Journal of Political Science 50.


# This file creates a figure that has one panel for each of the six outcomes.  
# Each panel plots four points -- one point for each model that I estimate.  


if (interactive()) {
  STANDARDIZE <- FALSE  # standardize the outcomes?
} else {
  clArgs <- commandArgs(trailingOnly = TRUE)
  if (length(clArgs) >= 1) { 
    STANDARDIZE <- as.logical(clArgs[1])  
  } else {
    writeLines(paste0("clArgs has length ", length(clArgs), ".  This length is less than 1, so the script is stopping."))
    stop()
  } 
}
  
library(Bullock, lib.loc = c(.libPaths(), 'packageLibrary'))  # for qw()
library(grid)
library(ivpack)   # for cluster.robust.se()
library(lattice)
library(Hmisc)    # for Dotplot()

dirOutput <- 'float_output/'  
source('IV_setup.R')
source('functions/estimateModels.R')
source('float_code/dotplotParameters.R')

filenameStem <- 'Figure_2'
PDFtitle     <- 'Figure 2: Effects of education on economic attitudes'
if (STANDARDIZE) { 
  filenameStem <- 'Figure_A05'
  PDFtitle     <- paste(PDFtitle, '(standardized outcomes)', PDFtitle) 
}
varNames    <- qw('eqwlth goveqinc guarantee.7pt govt.health.7pt helppoor welfare')
dfNames     <- qw('GSS.df GSS.df   ANES.df       ANES.df         GSS.df   GSS.df') 



##############################################################################
# PRELIMINARIES FOR DRAWING THE FIGURE
##############################################################################
panelHeight      <- list(1.10, "inches") 
yBetween         <- 2.3000  # space between rows
xBetween         <- 0.8500  # space between columns
panelWidth       <- list(1.14, "inches")        
panelBorderWidth <- .3
panelLayout      <- c(3, 2)  # columns, then rows
stripHeight      <- list(lines = 2.55, lineheight = 14) # lineheight is space between lines
nPanels          <- panelLayout[1] * panelLayout[2] 
baseCexSize      <- 1.00     # cex
stripTextSize    <-  .700 * baseCexSize 
xLabSize         <- 1.000 * baseCexSize 
yLabSize         <- 1.000 * baseCexSize
xAxisTextSize    <-  .670 * baseCexSize 
yAxisTextSize    <-  .710 * baseCexSize
axisTickSize     <- .4
if (STANDARDIZE) {
  xListLimits <- c(-.3, 1.2)        
  xListLabels <- qw("-.15 .25 .65 1.05")  
} else {    
  xListLimits <- c(-.15, .35)  
  xListLabels <- qw("-.1 0 .1 .2 .3")  
}
xListAt <- as.numeric(xListLabels)
xList <- list(
  draw        = TRUE,         
  limits      = xListLimits,
  labels      = xListLabels, 
  at          = xListAt,     
  tck         = c(axisTickSize, 0), 
  col         = panelBorderCol, 
  cex         = xAxisTextSize,
  alternating = 1, 
  relation    = 'free', 
  axs         = 'i')  # axs='i' means that there is no padding around the xlimits
yList <- list(
  draw        = TRUE,
  limits      = c(.5, 4.6),
  labels      = c(
    "Baseline model (Eq. 2)",
    "\U2026with cohort\U00ADyear fixed effects",
    "\U2026with state\U00ADwhen\U00ADyoung \U00D7 year\U00ADwhen\U00ADyoung trends",
    "\U2026with political controls",
    "\U2026with demographic controls",
    "\U2026with political and demographic controls"),
  at          = 4:1,
  tck         = 0,
  col         = panelBorderCol,
  cex         = yAxisTextSize,
  alternating = c(1, 1),
  relation    = 'same',  
  axs         = 'i')
if (xList$labels[1] == '-.1') {
  xList$at[1] <- -.115  # visual centering for a label w/ a negative sign
}
yList$labels <- yList$labels[c(1:3, 6)]



##############################################################################
# STANDARDIZE OUTCOMES 
##############################################################################
if (STANDARDIZE) {
  for (i in 1:length(varNames)) {
    myDF <- get(dfNames[i])
    myDF[, varNames[i]] <- scale(myDF[, varNames[i]])
    assign(dfNames[i], myDF)
  }
  rm(i, myDF)
}



##############################################################################
# ESTIMATE MODELS
##############################################################################
IVModelObjectsEnv <- estimateModels(
  varNames     = varNames, 
  modelEnvir   = as.environment(IVModelsEnv), 
  dfNames      = dfNames,
  objectSuffix = '.IV') 
IVClusters <- lapply(
  IVModelObjectsEnv,
  function (x) { 
    modelRHS <- as.character(x$formula[3]) 
    if (grepl('factor(yearYoungNorm)', modelRHS, fixed = TRUE)) { 
      droplevels(x$model$stateYoung:x$model$'factor(yearYoungNorm)')
    } else {
      droplevels(x$model$stateYoung:factor(x$model$yearYoungNorm))
    }
  }
)



##############################################################################
# FIND MINIMUM AND MAXIMUM SAMPLE SIZES
##############################################################################
nMin <- apply(
  X      = getNobsForMatrixOfModels(varNames, objectEnvir = IVModelObjectsEnv), 
  MARGIN = 2, 
  FUN    = min)
nMax <- apply(
  X      = getNobsForMatrixOfModels(varNames, objectEnvir = IVModelObjectsEnv), 
  MARGIN = 2, 
  FUN    = max)
nMin
nMax



##############################################################################
# CREATE DATA FRAME WITH DATA TO PLOT
##############################################################################
treatmentName <- 'educ'
dataToPlot <- expand.grid(
  lower     = NA, 
  pred      = NA, 
  upper     = NA,
  model     = 1:length(IVModelsEnv),
  depVar    = varNames)
for (varName in varNames) {
  varNameIndex <- which(varNames == varName)
  for (modNum in unique(dataToPlot$model)) {
    tmpModObj  <- get(paste0(varName, '.IV', modNum), IVModelObjectsEnv)
    tmpCluster <- get(paste0(varName, '.IV', modNum), IVClusters)
    coefEst    <- coeftest(tmpModObj)[treatmentName, 'Estimate']
    coefSE     <- cluster.robust.se(tmpModObj, tmpCluster)[treatmentName, 'Std. Error']
    dataToPlot[dataToPlot$model == modNum & dataToPlot$depVar == varName, c("lower", "pred", "upper")] <- c(
      coefEst - 1.96*coefSE,
      coefEst,
      coefEst + 1.96*coefSE)
  } 
}  

# Order the panels according to the order of varNames, which is the order that
# I use for the tables in the paper.  
for (i in 1:length(varNames)) {
  dataToPlot$panelNum[grepl(varNames[i], dataToPlot$depVar)] <- i
}

# Add row numbers, i.e., y-coordinates for each point that I plot.  For 
# example, if results from four models are to be plotted in each panel, 
# length(IVModelsEnv) == 4, and the assigned row numbers will be 1 through 4.
dataToPlot$rowNum <- length(IVModelsEnv)+1 - dataToPlot$model



##############################################################################
# STRIP FUNCTION
##############################################################################
horizStrip <- function(...) {
  stripYCoords    <- c(0,0,1,1,0)
  stripTextYCoord <- .5 
  lpolygon(
    x        = c(0,1,1,0,0), 
    y        = stripYCoords, 
    col      = stripBackgroundCol, 
    border   = stripBorderCol,
    linejoin = 'mitre',  # for square corners
    lwd      = panelBorderWidth)
  ltext(
    x          = .5, 
    y          = stripTextYCoord, 
    labels     = stripText[panel.number()], 
    font       = 2, 
    cex        = stripTextSize, 
    lineheight = 1)  # .4 for one-line labels         
}  



##############################################################################
# PANEL FUNCTION
##############################################################################
myPanel <- function(...) {
  panel.abline(v = 0, col = dotplotLineColor)
  panel.Dotplot(...)
  
  # Put tick marks at top of each panel
  panel.axis(
    side        = "top",
    at          = xList$at,
    draw.labels = FALSE,
    half        = FALSE,
    tck         = axisTickSize * .75)
}



##############################################################################
# MAKE THE PDF FILE
##############################################################################
pdf(
  file   = paste0(dirOutput, filenameStem, '.pdf'), 
  width  = PS_width, 
  height = PS_height, 
  paper  = "special", 
  title  = PDFtitle,
  bg     = postscriptBackground)
trellis.par.set(
  axis.components = list(     # distance between axis ticks and tick labels
    left   = list(pad1 = .50),
    bottom = list(pad1 = .80)),   
  axis.line = list(
    alpha = 1, 
    col   = panelBorderCol, 
    lty   = 1, 
    lwd   = panelBorderWidth),  
  dot.line = list(alpha=1, col=dotplotLineColor, lty=1, lwd=1),
  dot.symbol = list(
    alpha = 1, 
    cex   = .5, 
    col   = dotColor, 
    font  = 1, 
    pch   = 16),           
  plot.line = list(alpha = 1, col = dotColor, lty = 1, lwd = 1),
  superpose.line = list(
    alpha = c(1,1), 
    col   = c('black', plotLineColor), 
    lty   = c("solid", "43"), 
    lwd   = c(1, plotLineWidth)))
mainResultsDotplot <- Dotplot(
  rowNum ~ Cbind(pred, lower, upper) | panelNum, 
  data           = dataToPlot,
  panel          = myPanel,
  layout         = panelLayout,
  as.table       = TRUE,
  between        = list(x = xBetween, y = yBetween),                   
  strip          = horizStrip,
  par.strip.text = stripHeight,
  scales         = list(x = xList, y = yList),
  xlab           = '',
  ylab           = '',
  par.settings   = list(clip = list(panel="on", strip="on") ))
print(
  mainResultsDotplot, 
  panel.width  = panelWidth, 
  panel.height = panelHeight)


# SAVE AND PROCESS FILE
dev.off()
if (!(Sys.which('pdfcrop')=='' | Sys.which('perl')=='')  |  'pdfcrop' %in% installed.packages()[, 'Package']) {  # if "pdfcrop" is installed 
  system(paste(
    if (Sys.info()['sysname']=='Windows') paste(Sys.getenv('Comspec'), '/c ') else '',
    'pdfcrop', 
    paste0(dirOutput, filenameStem, '.pdf'), 
    paste0(dirOutput, filenameStem, '_crop.pdf')))
  file.remove(paste0(dirOutput, filenameStem, '.pdf'))
  if (interactive() & Sys.info()['sysname']=='Windows') shell.exec(paste0(normalizePath(dirOutput), paste0(filenameStem, '_crop.pdf')))
}